Self— Consistent Mode— Coupling Approach to ID Heat Transport 



O 

o 



o 

B 



X3 
O 

o 



> 

On 

m 
o 

o 



I 

o 
o 



X 
S3 



Luca Delfini/ Stefano Lepri,^'E] Roberto Livi,^'0 and Antonio Politi^ 

'istituto Nazionale di Ottica Applicata, largo E. Fermi 6 1-50125 Firenze, Italy 
^Istituto dei Sistemi Complessi, Consiglio Nazionale delle Ricerche, largo E. Fermi 6 1-50125 Firenze, Italy 
Dipartimento di Fisica, via G. Sansone 1 1-50019, Sesto Fiorentino, Italy 
(Dated: February 6, 2008) 

In the present Letter we present an analytical and numerical solution of the self-consistent mode- 
coupling equations for the problem of heat conductivity in one-dimensional systems. Such a solution 
leads us to propose a different scenario to accomodate the known results obtained so far for this 
problem. More precisely, we conjecture that the universality class is determined by the leading order 
of the nonlinear interaction potential. Moreover, our analysis allows us determining the memory 
kernel, whose expression puts on a more firm basis the previously conjectured connection between 
anomalous heat conductivity and anomalous diffusion. 



PACS numbers: 63.10.-|-a 05.60.-k 44.10.+i 

It is well known that relaxation and transport phe- 
nomena in reduced spatial dimensions {d < 3) are of- 
ten qualitatively different from their three-dimensional 
counterparts. This is a documented effect, for example, 
in single-filing systems, where particle diffusion does not 
follow Pick's law . Another related phenomenon is the 
enhancement of vibrational energy transmission in quasi- 
ID systems like polymers |3| or individual carbon nan- 
otubes 3]. The specific instance of anomalous thermal 
conduction in low-dimensional many-particle systems has 
recently received a renewed attention (3| . Anomalous be- 
haviour means both a divergence of the finite-size conduc- 
tivity k{L) (X L" in the large-size limit L — > oo and a 
nonintegrable decay of equilibrium correlations of the en- 
ergy current (the Green-Kubo integrand), (J(t)J(O)) oc 
^-(i-a) large times i ^ oo (with < a < 1). Sim- 
ulations jif and theoretical arguments [gj indicate that 
anomalies should occur generically in d < 2 whenever 
momentum is conserved. 

The importance of predicting the scaling behaviour 
(i.e. the value of a) is twofold: (i) on a basic ground, to 
classify the ingredients (e.g. symmetries) that define the 
possible universality classes; (ii) on a practical ground, 
to estimate heat conductivity in finite systems, a crucial 
issue to compare with experimental data on, say, carbon 
nanotubes. In spite of several efforts, the theoretical sce- 
nario is still controversial. In d = 1, arguments based 
on Mode-Coupling Theory (MCT) [1 a well-known 
approach to estimate long-time tails of fluids and to 
describe the glass transition yield a = 2/5. This 
estimate was criticized as inconsistent in Ref. where 
renormalization group arguments were instead shown to 
give a = 1/3. Nevertheless, the 2/5 value has been later 
derived both from a kinetic-theory calculation for the 
quartic (/3) Fermi Pasta Ulam (FPU) model [ll| and from 
a solution of the MCT by means of an ad hoc Ansatz 
[l^ . It was thereby conjectured that 2/5 is found 
for a purely longitudinal dynamics, while a crossover to- 
wards 1/3 can to be observed only in the presence of a 



coupling to transversal motion. Unfortunately, the accu- 
racy of numerical simulations is generally insufficient to 
disentangle the whole picture. The only two convincing 
studies concern the hard point gas, which has been re- 
cently found to be characterized by a — 1/3 ^'^d 
the purely quartic FPU model, where instead a is defi- 
nitely larger than 1/3 (and possibly closer to 2/5) [l3 |. 
The situation is even more controversial in d = 2 where 
logarithmic divergence is expected jl5j . 

The exact self-consistent solution of the MCT equa- 
tions presented in this Letter demonstrates that the over- 
all scenario is different from that one proposed in 1^, 
namely that a — 1/3 in the presence of cubic nonlin- 
earities. This prediction is confirmed by our numerical 
simulations of the FPU model with cubic potential which 
yields sizably different a values with respect to the quar- 
tic case. Altogether, theoretical and numerical results 
indicate that the asymptotic scaling behaviour is deter- 
mined by the order of the leading nonlinearity in the 
interaction potential. 

Let us consider the simplest one-dimensional version 
of the self-consistent MCT for the normalized correla- 
tor of the Fourier transform of the displacement field 
G{q,t) = {Q*{q,t)Q{q,0))/{\Q{q)\^). In dimensionless 
units in which the particle mass, the lattice spacing and 
the bare sound velocity are set to unity, they read [a, 0| 



Giq,t) 

T{q,t)^ u;-\q) 



e I T{q, t - s)G{q, s) ds -f Lj'iq)Giq, t) = 


^^-^^'^ Gip,t)G(p',t) 
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p+p' — g— 0,±7r 



We consider periodic boundaries so that the wavenum- 
bers are q = 27rfc/7V with -N/2 +l<k< N/2. Notice 
that G{q, t) = G{—q, t). Eqs. must be solved with the 
initial conditions G{q, 0) = 1 and G{q, 0) — 0. 

The first of Eqs. Q is exact and is derived within the 
well-known Mori-Zwanzig projection approach In 
the small-wavenumber limit, it describes the response of 
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an elastic string at finite temperature. The above mode- 
coupling approximation of the memory function T has 
been derived for a chain of atoms interacting through 
a nearest-neighbour anharmonic potential V{x) 0, 0| 
whose expansion around its minimum at a: = is of the 
form x"^ /2+g3x^ /3 + . . . (see, e.g., the Lennard- Jones po- 
tential). Both the coupling constant s and the dispersion 
relation io{q) are temperature-dependent input parame- 
ters that must be computed independently by simulation 
or approximate analytical approaches [Sl, ll^l- For the 
aims of the present Letter, we may restrict ourselves to 
considering their bare values, obtained in the harmonic 
approximation that, in our units, read e = 3(/|fcBT/27r 
and u!{q) = 2|sin||. Of course, the actual renormal- 
ized values are needed when a quantitative comparison 
with a specific model is looked for. Moreover, since the 
anomalies we are interested in stem from the nonlinear 
interaction of long- wavelength modes, we let uj(q) — \q\, 
in the analytic treatment presented below. 

Direct numerical simulations ^ indicate that nonlin- 
ear and nonlocal losses in Eq. Q are small compared 
to the oscillatory terms. This suggests splitting the G 
dynamics into phase and amplitude evolution, 

G{q,t) = C(g,i)e'"^'^* +C.C. (2) 

Upon substituting this equation into Eq. Q, one obtains, 
in the slowly varying envelope approximation, qC 3> C, 

2^C((Z,t)+e^ dt'Miq,t-t')C{q,t') = (3) 

plus a similar expression for C* , while the new kernel M 
turns out to be 

oo 

Miq,t)^q^ J dpC*ip~q,t)Cip,t) , (4) 
— oo 

where the sum in has been replaced by a suitable 
integral, since we consider the thermodynamic limit N = 
oo and small q- values, which are, by the way, responsible 
for the asymptotic behavior. 

Notice that Eqs. (|3I4II have been obtained after dis- 
carding the second order time derivative ofC{q, t) as well 
as the integral term proportional to (7, besides all rapidly 
rotating terms. The validity of this approximation is re- 
lated to the separation between the decay rate of C(g, t) 
and uj{q)\ its correctness will be cheked a posteriori, after 
discussing the scaling behaviour of C{q,t). Notice also 
that in this approximation, Umklapp processes do not 
contribute: it is in fact well known that they are negligi- 
ble for long-wavelength phonons in ID Q . 

Having transformed the second order differential equa- 
tion for G into a first order one for C, we can introduce a 
simple scaling argument yielding the dependence of G on 
q and t as follows (see also , where a similar equation 



was investigated), 

Giq,t) = giV^tq'/') , M{q,t) = q^ f{V^tq^/^). (5) 

This shows that the decay rate for the evolution oiG{q, t) 
is of the order q'^^'^^/e, which has to be compared with the 
scale q of the corresponding phase factor. Accordingly, 
amplitude and phase dynamics become increasingly sep- 
arated for q {). High q- values (g w 1) are those for 
which the slowly varying envelope approximation is less 
accurate. However, if e is small enough, such modes are 
correctly described, too. 

The functions / and g can be determined by substitut- 
ing the previous expressions into Eqs. H3I4|I . Upon setting 
X — \fetq^^'^ , one obtains the equation 

g'(x) = -\ dyf{x-y)g{y) (6) 
Jo 

/ + 00 
dyg*{\xy^-y\y')g{yy'){7) 
-oo 

The asymptotic behaviour for x —^ can be determined 
analytically, 

, , 1 / ax^/'^x , a 

5(x) = -exp(^~^j , f{x) ^ (8) 

where a is a suitable constant which is determined self- 
consistently from Eq. iQ. To assess the validity of 
the above calculation, we have numerically integrated 
Eqs. iQJ by the Euler method for the original dispersion 
relation u}{q) and different e- values. We have verified 
that a time step At = 0.01 guarantees a good numeri- 
cal accuracy over the explored time range. The Fourier 
transform G{q, uj) is plotted in Fig. ^ for three different 
q- values versus w — LL>maxiq), where ujmax{q) is the fre- 
quency corresponding to the maximal value Gmax of the 
spectrum (this is equivalent to removing the oscillating 
component from G{q,t)). Furthermore, in order to test 
relation (jSJ, the vertical axis is scaled to the maximum G- 
value, while the frequencies are divided by the half-width 
j{q) at half of the maximum height. This latter quantity 
can be interpreted as the inverse lifetime of fluctuations 
of wavenumber q. The good data collapse confirms the 
existence of a scaling regime. The approximate analyti- 
cal expression is in excellent agreement with the data. As 
expected, some deviations are present for small lo where 
Eq. (|SJ) is not strictly applicable. Moreover, in the in- 
set of Fig. n where the same curves are plotted using 
doubly logarithmic scales, one sees that the lineshapes 
follow the predicted power-law, w"''/'^, over a wide range 
of frequencies. In Fig. |21 we show that j{q) is indeed 
proportional to ^/e (f"!^ . It is particularly instructive to 
notice that the agreement is very good also also for a 
relatively large value of the coupling constant (e ~ 1), 
although the slowly varying envelope approximation is 
not correct for large q values. The deviations observed 
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FIG. 1: Fourier transform G(q, uj) of the correlation functions 
for 3 different wavenumbers e = 1, N = 2000. The solid line is 
the lineshape computed by the approximate analytic theory, 
i.e. by Fourier transforming the function C(g, t) defined in 
Eqs. I5I8II . The same curves are plotted in the inset in log-log 
scales, where only positive frequencies are shown. 



FIG. 3: The logarithmic derivative of the energy-flux spec- 
trum ^(i^) versus the frequency v for the FPU potential with 
energy density set equal to 10. The two horizontal lines cor- 
respond to the theoretical predictions —1/3 and —2/5. The 
statistical error is on the order of the observed irregular fluc- 
tuations. 
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FIG. 2: Scaling of the linewidth 7(g) of G{q,Lu) with q for 
3 different values of the coupling constant e and A'' = 2000. 
The solid line corresponding to the power law cf'^^ is plotted 
for reference. 



at small g-values for small couplings are due to the very 
slow convergence in time. Better performances could be 
obtained by increasing both N then the integration time 
(10^, in our units) well beyond our current capabilities. 

It is crucial to compare these result with the previ- 
ous work. Making use of Eqs. H5I8|I . it can be shown 
that the memory function T contains terms of the form 
q2 g±jgt^^2/3^ j j-j- oscillates with a power-law envelope. 
Accordingly, its Laplace trasform has branch-cut singu- 
larities of the form /{z ± q)3 . This finding is not con- 
sistent with the heuristic assumption of Rcfs. [3,|1| and 
the result of i where MCT equations were solved with 
the Ansatz V{q,z) = q^V{z). In addition, the numerical 
solution does not show any signature of the / z^/^ de- 
pendence found in |l2l| . For instance, it would imply a 
peak at w = in the spectrum of T which is, instead. 



absent in numerical solutions. 

In order to estimate the long-time decay of the energy- 
current autocorrelation we use the approximate expres- 
sion proposed in [l^ . 

{J{t)jm oc Y.v^{q)G\q,t) (9) 

where v{q) — w'{q) is the bare group velocity. Neglecting 
the oscillating terms in and letting v{q) ~ 1 (the sum 
is dominated by the small-q terms), we can use Eq. ^ 
to find 

(J(i)J(O)) cx J dqg^i^tq^/^) cc ^ (10) 

i.e. a = 1/3. Note that this result is independent of the 
actual form of g (provided convergence of the integral 
is insured). This scaling has also been checked to hold 
by directly evaluating the sum with the numerically 
computed correlations G. 

As already mentioned at the beginning, several numerical 
simulations have been performed to estimate the expo- 
nent a. Since the most accurate data Q| differ signif- 
icantly from 1/3, we decided to run a new set of sim- 
ulations. We considered the FPU model with interpar- 
ticle potential V{x) = + gsx^/S + x'^/A and peri- 
odic boundary conditions. We performed microcanoni- 
cal simulations to cornpute the average power spectrum 
S{iy) of J (see Refs. [1 E3 for details). The long-time 
tail manifests itself as a power-law divergence at 
low- frequencies, v 0. In order to provide a reliable es- 
timate of a, it is convenient to evaluate the logarithmic 
derivative d\n S/dlni^. As shown in Fig. 01 for 33 = 1 this 
quantity does display a plateau around —1/3 (the growth 
towards zero at very small v values is due to the cutoff 
introduced by the finite size of the lattice) . On the other 
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hand, the data obtained for 53 = (FPU-/3 model) indi- 
cate a sensibly different scaling exponent, which is much 
closer to —2/5 and in agreement with previous works 

We have thus reached the important conclusion that 
the memory kernel decays algebraically in the prescribed 
regime and, accordingly, the relaxation is not exponen- 
tial logC ~ -^2^4/3 (i non Lorenzian lineshapes). The 
further striking feature is that conventional hydrodynam- 
ics breaks down, since the peak widths scale as g^^^ rather 
than q^, as expected in the standard case. In addition, 
the linewidths are connected to transport coefficients be- 
ing proportional to Ag^, where A is the sound attenua- 
tion constant. The anomalous scaling can be recasted in 
terms of a diverging A{q) ^ q^^^^. Altogether, one may 
think of this as a superdiffusive process, intermediate be- 
tween standard diffusive and ballistic propagation. Our 
result thus strenghten the picture emerging in Ref. 
from the analysis of the hard-point gas, where it has 
been shown that energy perturbations perform a Levy 
walk. One merit of our approach is that it allows for 
a direct connection with anomalous diffusion problems 
[l9| . As it is known, these can be modeled by generalized 
Langevin equations with power-law kernels. If we now as- 
sume that expression (jHl for / holds for every x, we can 
solve Eqs. (jSJ by Laplace transforming Eq. to obtain 
C{q,z) = iz^l'^ I {iz^l''^ -|- aq^\ This expression is pre- 
cisely the Laplace transform of the Mittag-Leffler func- 
tion Ef,[-{\tY) Uli^ for /z = 4/3 and A = {aq^fl^ 
|2ll |. This observation suggests that the effective evolu- 
tion of fluctuations should be modeled by the fractional 
differential equation 

|5^C((7,i) + A^C(<z,0 = 0. (11) 

The case of interest here (1 < /i < 2) corresponds to the 
so-called fractional oscillations ]20j. It should be em- 
phasized that in the present context, memory arises as a 
genuine many-body effect and needs not to be postulated 
a 'priori. 

In conclusion, we have shown that MCT with cubic 
nonlinearity predicts a decay of the heat current 
autocorrelation, i.e. a — 1/3. Our analysis reconciles 
this approach with the renormalization-group calculation 
Is] and supports the idea that the mechanisms yielding 
anomalous transport in ID are largely universal. The 
sizeable deviations observed for quartic potentials sug- 
gest the existence of a different universality class which 
should be described by different mode-coupling equations 
with a quartic nonlinearity. A preliminary analysis con- 
firms that this scenario is indeed correct 22] . Finally, we 
have analytically shown that memory effects emerging 
from the nonlinear interaction of long-wavelength modes 
can be described by a generalized Langevin equation with 
power-law memory. This provides a sound basis estab- 



lishing a connection between anomalous transport and 
superdiffusive processes. 
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